#pragma once


#include <iostream>
#include <list>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <filesystem>
#include <unordered_map>
#include <algorithm>
#include <functional>
#include <cstdlib>
#include <cmath>


#include <string>
#include <vtkExodusIIWriter.h>
#include <vtkSmartPointer.h>
#include <vtkPolyData.h>


#include <vtkType.h>
#include <vtkPoints.h>
#include <vtkPointData.h>
#include <vtkFloatArray.h>
#include <vtkHexahedron.h>
#include <vtkCellArray.h>
#include <vtkExplicitStructuredGrid.h>
#include <vtkDataSetMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkScalarBarWidget.h>
#include <vtkScalarBarActor.h>
#include <vtkLookupTable.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkAutoInit.h>
#include <vtkTetra.h>
#include <vtkTriangle.h>
#include <vtkPolygon.h>
#include <vtkTriangleStrip.h>

#include <vtkType.h>
#include <vtkPoints.h>
#include <vtkPointData.h>
#include <vtkFloatArray.h>
#include <vtkHexahedron.h>
#include <vtkCellArray.h>
#include <vtkExplicitStructuredGrid.h>
#include <vtkDataSetMapper.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkScalarBarWidget.h>
#include <vtkScalarBarActor.h>
#include <vtkLookupTable.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkAutoInit.h>
#include <vtkTetra.h>
#include <vtkTriangle.h>
#include <vtkPolygon.h>
#include <vtkTriangleStrip.h>
#include <vtkXMLUnstructuredGridWriter.h>
#include <vtkExodusIIWriter.h>

#include <vtkTimeSourceExample.h>
//#include "exodusII.h"
#include <inttypes.h>
#include <vtkExodusIIWriter.h>

#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkCompositeDataGeometryFilter.h>
#include <vtkDataSetMapper.h>
#include <vtkExodusIIReader.h>
#include <vtkGaussianKernel.h>
#include <vtkInterpolateDataSetAttributes.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkNamedColors.h>
#include <vtkPointData.h>
#include <vtkPointInterpolator.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkStaticPointLocator.h>
#include <vtkUnstructuredGrid.h>

#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkContourFilter.h>
#include <vtkImplicitModeller.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkVectorText.h>

#include <vtkImplicitModeller.h>

#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/features/normal_3d.h>
#include <pcl/surface/gp3.h>   
#include<pcl/io/vtk_lib_io.h>
#include <vtkNew.h>

#include <pcl/pcl_macros.h>
#include <pcl/PCLPointCloud2.h>
#include <pcl/PolygonMesh.h>

#include <pcl/io/vtk_lib_io.h>

#include <pcl/io/pcd_io.h>
#include <pcl/visualization/pcl_visualizer.h>

#include <vtkPolyDataWriter.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkUnstructuredGridWriter.h>
#include <vtkStructuredPointsWriter.h>

using namespace std;
struct GroupPoint3D {
   
    double xyz[3];
    int group_id;

    GroupPoint3D(double x, double y, double z, int group_id)
    {
        this->xyz[0] = x;
        this->xyz[1] = y;
        this->xyz[2] = z;
        this->group_id = group_id;
    }
    GroupPoint3D(double x, double y, double z)
    {
        this->xyz[0] = x;
        this->xyz[1] = y;
        this->xyz[2] = z;
    }
    void setValue(GroupPoint3D point)
    {
        xyz[0] = point.xyz[0];
        xyz[1] = point.xyz[1];
        xyz[2] = point.xyz[2];
        group_id = point.group_id;
    }
};
struct Interplator {

    /**
    vtkIdType vtkGaussianKernel::ComputeWeights(
        double x[3], vtkIdList* pIds, vtkDoubleArray* prob, vtkDoubleArray* weights)
    {
        vtkIdType numPts = pIds->GetNumberOfIds();
        double d2, y[3], sum = 0.0;
        weights->SetNumberOfTuples(numPts);
        double* p = (prob ? prob->GetPointer(0) : nullptr);
        double* w = weights->GetPointer(0);
        double f2 = this->F2;

        for (vtkIdType i = 0; i < numPts; ++i)
        {
            vtkIdType id = pIds->GetId(i);
            this->DataSet->GetPoint(id, y);
            d2 = vtkMath::Distance2BetweenPoints(x, y);

            if (vtkMathUtilities::FuzzyCompare(
                d2, 0.0, std::numeric_limits<double>::epsilon() * 256.0)) // precise hit on existing point
            {
                pIds->SetNumberOfIds(1);
                pIds->SetId(0, id);
                weights->SetNumberOfTuples(1);
                weights->SetValue(0, 1.0);
                return 1;
            }
            else
            {
                w[i] = (p ? p[i] * exp(-f2 * d2) : exp(-f2 * d2));
                sum += w[i];
            }
        } // over all points

        // Normalize
        if (this->NormalizeWeights && sum != 0.0)
        {
            for (vtkIdType i = 0; i < numPts; ++i)
            {
                w[i] /= sum;
            }
        }

        return numPts;
    }

    vtkIdType vtkProbabilisticVoronoiKernel::ComputeWeights(
        double x[3], vtkIdList* pIds, vtkDoubleArray* prob, vtkDoubleArray* weights)
    {
        vtkIdType numPts = pIds->GetNumberOfIds();
        double* p = (prob ? prob->GetPointer(0) : nullptr);
        double highestProbability = VTK_FLOAT_MIN;
        vtkIdType id, mostProbableId = 0;

        if (p) // return the point in the neighborhood with the highest probability
        {
            for (vtkIdType i = 0; i < numPts; ++i)
            {
                if (p[i] > highestProbability)
                {
                    mostProbableId = pIds->GetId(i);
                    highestProbability = p[i];
                }
            }
        }

        else // return the closest point in the footprint provided
        {
            double y[3], d, minD = VTK_FLOAT_MAX;
            for (vtkIdType i = 0; i < numPts; ++i)
            {
                id = pIds->GetId(i);
                this->DataSet->GetPoint(id, y);
                d = vtkMath::Distance2BetweenPoints(x, y);
                if (vtkMathUtilities::FuzzyCompare(d, 0.0,
                    std::numeric_limits<double>::epsilon() * 256.0)) // precise hit on existing point
                {
                    mostProbableId = id;
                    break;
                }
                else if (d <= minD)
                {
                    mostProbableId = id;
                    minD = d;
                }
            } // over all points
        }

        // Okay let's get out
        pIds->SetNumberOfIds(1);
        pIds->SetId(0, mostProbableId);
        weights->SetNumberOfTuples(1);
        weights->SetValue(0, 1.0);

        return 1;
    }
    */
    /*
    void writeExodusII(string filename) {

        vtkNew<vtkExodusIIWriter> vtkExodusIIWriter ;
        vtkExodusIIWriter->SetFileName(filename.c_str());
        vtkExodusIIWriter->WriteOutBlockIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalNodeIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalElementIdArrayOn();
        vtkExodusIIWriter->WriteAllTimeStepsOn();
        //print "#time", vtkTemporalDataSet.GetNumberOfTimeSteps();
        vtkExodusIIWriter->SetInputData(vtkTemporalDataSet);
        vtkExodusIIWriter->Update();
        vtkExodusIIWriter->Write();
    }*/
    
    void writeExodusII_A()
    {
        vtkNew<vtkTimeSourceExample> timeSource;

        vtkNew<vtkExodusIIWriter> exodusWriter;
        exodusWriter->SetFileName("output.e");
        exodusWriter->SetInputConnection(timeSource->GetOutputPort());
        exodusWriter->WriteAllTimeStepsOn();
        exodusWriter->Write();
    }

    void writeExodus()
    {
        // vtkNew<vtkTimeSourceExample> timeSource

        float x[8][3] =
        {
            {0,0,0},{1,0,0},{1,1,0},{0,1,0},
            {0,0,1},{1,0,1},{1,1,1},{0,1,1}
        };

        vtkNew<vtkPolyData> geometry;
        vtkNew<vtkPoints> points;
        vtkNew<vtkCellArray> polys;

        polys->InsertNextCell(8);
        for (size_t i = 0; i < 8; i++)
        {
            points->InsertPoint(i, x[i]);
            polys->InsertCellPoint(i);
        }

        geometry->SetPoints(points);
        geometry->SetVerts(polys);

        ;
        vtkNew<vtkPolyData> polyDataSource;
        vtkNew<vtkExodusIIWriter> exodusWriter;
        exodusWriter->SetFileName("result/polydata.e");
        //exodusWriter->SetInputConnection(geometry->GetOutput());
        //exodusWriter->SetInputData(geometry);
        exodusWriter->WriteAllTimeStepsOn();
        exodusWriter->Write();
    }

    void poly_triangle()
    {
        float x[8][3] =
        {
            {0,0,0},{1,0,0},{1,1,0},{0,1,0},
            {0,0,1},{1,0,1},{1,1,1},{0,1,1}
        };
        vtkIdType pts[6][4] =
        {
            {0,1,2,3},{4,5,6,7},{0,1,5,4},
            {1,2,6,5},{2,3,7,6},{3,0,4,7}
        };
        vtkNew<vtkPolyData> geometry;
        vtkNew<vtkPoints> points;
        for (size_t i = 0; i < 8; i++)
        {
            points->InsertPoint(i, x[i]);

        }
        vtkNew<vtkCellArray> polys;
        for (size_t i = 0; i < 6; i++)
        {
            polys->InsertNextCell(3, pts[i]);

        }

        geometry->SetPoints(points);
        //geometry->SetPolys(polys);
        geometry->SetStrips(polys);
        //geometry->SetVerts(polys);

        vtkNew<vtkPolyDataMapper> geometryMapper;
        geometryMapper->SetInputData(geometry);
        vtkNew<vtkActor> geometryActor;
        geometryActor->SetMapper(geometryMapper);

        vtkNew<vtkRenderer> renderer;
        renderer->AddActor(geometryActor);
        renderer->ResetCamera();
        renderer->SetBackground(0, 0, 0);

        vtkNew<vtkRenderWindow> renWin;
        renWin->AddRenderer(renderer);
        renWin->SetSize(300, 300);

        vtkNew<vtkRenderWindowInteractor> iren;
        iren->SetRenderWindow(renWin);

        renWin->Render();
        iren->Start();

    }

    vtkSmartPointer<vtkPolyData>  readYoungProp(const string& file_name) {

        fstream fin;

        fin.open(file_name, ios::in);

        if (!fin.is_open())
        {
            cout << "File not found !" << file_name << endl;
        }
        string line;
        double x, y, z;
        vector<GroupPoint3D> points_v;

        while (getline(fin, line)) {
            int i = 0;
            string arr[3];
            stringstream ssin(line);
            string arr_;
            while (ssin.good() && i < 3) {
                ssin >> arr[i];
                if (i == 0) {
                    arr_ = arr[0].substr(2);
                    x = strtod(arr_.c_str(), NULL);
                }
                else if (i == 1) {
                    arr_ = arr[1].substr(2);
                    y = strtod(arr_.c_str(), NULL);
                }
                else if (i == 2) {
                    arr_ = arr[2].substr(2);
                    z = strtod(arr_.c_str(), NULL);
                }
                i++;
            }
            GroupPoint3D gp3d = GroupPoint3D(x, y, z);
            points_v.push_back(gp3d);
        }
        cout << "size :" << points_v.size() << endl;
        fin.close();

        /////////////////////////////////////////////////////////
        vtkSmartPointer<vtkNamedColors> colors = vtkSmartPointer<vtkNamedColors>::New();

        // vtkNew<vtkUnstructuredGrid> ugrid = vtkNew<vtkUnstructuredGrid>::vtkNew();
        vtkNew<vtkUnstructuredGrid> ugrid;
        // Setup four points
        vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
        int i = 0;
        for (auto a : points_v) {
            //points->InsertNextPoint(a.xyz);
            points->InsertPoint(i, a.xyz);
            i++;
        }
        cout << "size : " << points_v.size() << endl;

        // Create the polygon
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        //vtkNew< vtkPolygon> polygon;
        polygon->GetPointIds()->SetNumberOfIds(18080); //make a quad

        double pts[4520][4];

        int in = 0;
        for (int i = 0; i < 4520; i++) {
            for (int s = 0; s < 4; s++) {
                pts[i][s] = in + s;
            }
            in = in + 4;
        }
        for (int i = 0; i < 4520; i++) {
            for (int j = 0; j < 4; j++) {
                polygon->GetPointIds()->SetId(pts[i][j], i);
            }
        }
        for (int i = 0; i < 4520; i++) {
            vtkSmartPointer <vtkIdList > idList_v = vtkSmartPointer<vtkIdList>::New();
            for (int j = 0; j < 4; j++) {
                idList_v.Get()->InsertNextId(pts[i][j]);
            }
            ugrid->InsertNextCell(VTK_QUAD, idList_v);
        }
        ugrid->SetPoints(points);

        // Add the polygon to a list of polygons
        vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();
        polygons->InsertNextCell(polygon);

        // Create a PolyData
        vtkSmartPointer<vtkPolyData> polygonPolyData = vtkSmartPointer<vtkPolyData>::New();
        polygonPolyData->SetPoints(points);
        polygonPolyData->SetPolys(polygons);
        polygonPolyData->BuildLinks();

        
        return polygonPolyData;
    }

    void buildYoungmodQuadWithhUnstructuredGrid(const string& file_name, const string& output_name)
    {
        fstream fin;
        fin.open(file_name, ios::in);

        if (!fin.is_open())
        {
            cout << "File not found !" << file_name << endl;
        }
        string line;
        double x, y, z;
        vector<GroupPoint3D> points_v;

        while (getline(fin, line)) {
            int i = 0;
            string arr[3];
            stringstream ssin(line);
            string arr_;
            while (ssin.good() && i < 3) {
                ssin >> arr[i];
                if (i == 0) {
                    arr_ = arr[0].substr(2);
                    x = strtod(arr_.c_str(), NULL);
                }
                else if (i == 1) {
                    arr_ = arr[1].substr(2);
                    y = strtod(arr_.c_str(), NULL);
                }
                else if (i == 2) {
                    arr_ = arr[2].substr(2);
                    z = strtod(arr_.c_str(), NULL);
                }
                i++;
            }
            GroupPoint3D gp3d = GroupPoint3D(x, y, z);
            points_v.push_back(gp3d);
        }
        cout << "size :" << points_v.size() << endl;
        fin.close();

        /////////////////////////////////////////////////////////
        vtkSmartPointer<vtkNamedColors> colors = vtkSmartPointer<vtkNamedColors>::New();

        // vtkNew<vtkUnstructuredGrid> ugrid = vtkNew<vtkUnstructuredGrid>::vtkNew();
        vtkNew<vtkUnstructuredGrid> ugrid;
        // Setup four points
        vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
        int i = 0;
        for (auto a : points_v) {
            //points->InsertNextPoint(a.xyz);
            points->InsertPoint(i, a.xyz);
            i++;
        }
        cout << "size : " << points_v.size() << endl;

        // Create the polygon
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        //vtkNew< vtkPolygon> polygon;
        polygon->GetPointIds()->SetNumberOfIds(18080); //make a quad

        double pts[4520][4];

        int in = 0;
        for (int i = 0; i < 4520; i++) {
            for (int s = 0; s < 4; s++) {
                pts[i][s] = in + s;
            }
            in = in + 4;
        }
        for (int i = 0; i < 4520; i++) {
            for (int j = 0; j < 4; j++) {
                polygon->GetPointIds()->SetId(pts[i][j], i);
            }
        }
        for (int i = 0; i < 4520; i++) {
            vtkSmartPointer <vtkIdList > idList_v = vtkSmartPointer<vtkIdList>::New();
            for (int j = 0; j < 4; j++) {
                idList_v.Get()->InsertNextId(pts[i][j]);
            }
            ugrid->InsertNextCell(VTK_QUAD, idList_v);
        }
        ugrid->SetPoints(points);
        

        // Add the polygon to a list of polygons
        vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();
        polygons->InsertNextCell(polygon);

        // Create a PolyData
        vtkSmartPointer<vtkPolyData> polygonPolyData = vtkSmartPointer<vtkPolyData>::New();
        polygonPolyData->SetPoints(points);
        polygonPolyData->SetPolys(polygons);
        polygonPolyData->BuildLinks();
        // polygon->Delete();

        vtkNew<vtkDataSetMapper> ugridMapper;
        ugridMapper->SetInputData(ugrid);


        string exodusii_name = output_name;
        //vtkSmartPointer<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkNew<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkExodusIIWriter->SetFileName(output_name.c_str());
        vtkExodusIIWriter->WriteOutBlockIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalNodeIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalElementIdArrayOn();
        vtkExodusIIWriter->WriteAllTimeStepsOn();
        //print "#time", vtkTemporalDataSet.GetNumberOfTimeSteps();
        //vtkExodusIIWriter->SetInputData(polygonPolyData);
        vtkExodusIIWriter->SetInputData(ugrid);
        //vtkExodusIIWriter->SetInputConnection(ugridMapper->GetOutputPort());
        vtkExodusIIWriter->Update();
        vtkExodusIIWriter->Write();
        //vtkExodusIIWriter->Delete();

        vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
        actor->SetMapper(ugridMapper);
        actor->GetProperty()->SetColor(colors->GetColor3d("Silver").GetData());

        // Visualize
        vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
        vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
        renderWindow->SetWindowName("Polygon");
        renderWindow->AddRenderer(renderer);
        vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
        renderWindowInteractor->SetRenderWindow(renderWindow);

        renderer->AddActor(actor);
        renderer->SetBackground(colors->GetColor3d("Salmon").GetData());
        renderWindow->Render();
        renderWindowInteractor->Start();
    }


    void buildPissonQuadWithhUnstructuredGrid(const string& file_name, const string& output_name)
    {
        fstream fin;
        fin.open(file_name, ios::in);

        if (!fin.is_open())
        {
            cout << "File not found !" << file_name << endl;
        }
        string line;
        double x, y, z;
        vector<GroupPoint3D> points_v;

        while(getline(fin, line)){
            int i = 0;
            string arr[3];
            stringstream ssin(line);
            string arr_;
            while (ssin.good() && i < 3) {
                ssin >> arr[i];
                if (i == 0) {
                    arr_ = arr[0].substr(2);
                    x = strtod(arr_.c_str(), NULL);
                }
                else if (i == 1) {
                    arr_ = arr[1].substr(2);
                    y = strtod(arr_.c_str(), NULL);
                }
                else if (i == 2) {
                    arr_ = arr[2].substr(2);
                    z = strtod(arr_.c_str(), NULL);
                }
                i++;
            }
           GroupPoint3D gp3d = GroupPoint3D(x, y, z);
           points_v.push_back(gp3d);
        }
        cout << "size :" << points_v.size() << endl;
        fin.close();

       /////////////////////////////////////////////////////////
        vtkSmartPointer<vtkNamedColors> colors = vtkSmartPointer<vtkNamedColors>::New();

       // vtkNew<vtkUnstructuredGrid> ugrid = vtkNew<vtkUnstructuredGrid>::vtkNew();
        vtkNew<vtkUnstructuredGrid> ugrid;
        // Setup four points
        vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
        int i = 0;
        for (auto a : points_v) {
            //points->InsertNextPoint(a.xyz);
            points->InsertPoint(i, a.xyz);
            i++;
        }
        cout << "size : " << points_v.size() << endl;

        // Create the polygon
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        //vtkNew< vtkPolygon> polygon;
        polygon->GetPointIds()->SetNumberOfIds(18156); //make a quad

        double pts[4539][4];

        int in = 0;
        for (int i = 0; i < 4539; i++) {
            for (int s = 0; s < 4; s++) {
                pts[i][s] = in + s;
            }
            in = in + 4;
        }
        for (int i = 0; i < 4539; i++) {
            for (int j = 0; j < 4; j++) {
                polygon->GetPointIds()->SetId(pts[i][j], i);
            }
        }
        for (int i = 0; i < 4539; i++) {
            vtkSmartPointer <vtkIdList > idList_v = vtkSmartPointer<vtkIdList>::New();
            for (int j = 0; j < 4; j++) {
                idList_v.Get()->InsertNextId(pts[i][j]);
            }
            ugrid->InsertNextCell(VTK_QUAD, idList_v);
        }
        ugrid->SetPoints(points);

        
        // Add the polygon to a list of polygons
        vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();
        polygons->InsertNextCell(polygon);

        // Create a PolyData
        vtkSmartPointer<vtkPolyData> polygonPolyData = vtkSmartPointer<vtkPolyData>::New();
        polygonPolyData->SetPoints(points);
        polygonPolyData->SetPolys(polygons);
        polygonPolyData->BuildLinks();
        // polygon->Delete();

        vtkNew<vtkDataSetMapper> ugridMapper;
        ugridMapper->SetInputData(ugrid);

        string exodusii_name = output_name;
        //vtkSmartPointer<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkNew<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkExodusIIWriter->SetFileName(output_name.c_str());
        vtkExodusIIWriter->WriteOutBlockIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalNodeIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalElementIdArrayOn();
        vtkExodusIIWriter->WriteAllTimeStepsOn();
        //print "#time", vtkTemporalDataSet.GetNumberOfTimeSteps();
        //vtkExodusIIWriter->SetInputData(polygonPolyData);
        vtkExodusIIWriter->SetInputData(ugrid);
        //vtkExodusIIWriter->SetInputConnection(ugridMapper->GetOutputPort());
        vtkExodusIIWriter->Update();
        vtkExodusIIWriter->Write();
        //vtkExodusIIWriter->Delete();

        vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
        actor->SetMapper(ugridMapper);
        actor->GetProperty()->SetColor(colors->GetColor3d("Silver").GetData());

        // Visualize
        vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
        vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
        renderWindow->SetWindowName("Polygon");
        renderWindow->AddRenderer(renderer);
        vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
        renderWindowInteractor->SetRenderWindow(renderWindow);

        renderer->AddActor(actor);
        renderer->SetBackground(colors->GetColor3d("Salmon").GetData());
        renderWindow->Render();
        renderWindowInteractor->Start();
    }

    void buildQuadWithUnstructuredGrid(const string& file_name)
    {
        fstream fin;
        fin.open(file_name, ios::in);

        if (!fin.is_open())
        {
            cout << "File not found !" << file_name << endl;
        }
        string line;
        double x, y, z;
        int group_id;
        vector<GroupPoint3D> points_v;

        while (fin >> x >> y >> group_id >> z) {
           
            GroupPoint3D gp3d = GroupPoint3D(x, y, z, group_id);
            points_v.push_back(gp3d);
        }
        fin.close();

        ///////////////////////////////////////////////////////////////////////////////
        vtkSmartPointer<vtkNamedColors> colors = vtkSmartPointer<vtkNamedColors>::New();

        vtkNew<vtkUnstructuredGrid> ugrid = vtkNew<vtkUnstructuredGrid>::vtkNew();
        //vtkNew<vtkUnstructuredGrid> ugrid;
        // Setup four points
        vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
        int i = 0;
        for (auto a : points_v) {
            //cout << a.group_id << " " << x << " " << y << " " << z << endl;
            //points->InsertNextPoint(a.xyz);
            points->InsertPoint(i, a.xyz);
            i++;
        }

        // Create the polygon
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        //vtkNew< vtkPolygon> polygon;
        polygon->GetPointIds()->SetNumberOfIds(10236); //make a quad

        double pts[2559][4];
        
        int in = 0;
        for (int i = 0; i < 2559; i++) {
            for (int s = 0; s < 4; s++) {
                pts[i][s] = in+s;               
            }           
            in = in + 4;
         }
        for (int i = 0; i < 2559; i++) {
           for (int j = 0; j < 4; j++) {             
               polygon->GetPointIds()->SetId(pts[i][j], i);
           }
        }      
       for (int i = 0; i < 2559; i++) {
           vtkSmartPointer <vtkIdList > idList_v = vtkSmartPointer<vtkIdList>::New();
           for (int j = 0; j < 4; j++) {
               idList_v.Get()->InsertNextId(pts[i][j]);
           }
           ugrid->InsertNextCell(VTK_QUAD, idList_v);           
       }     
        ugrid->SetPoints(points);
        
        /*
        string fileName = "result/xmlunstructuredgrid.vtu";
        vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
        //writer.SetFileName('cosipy.vtu')
        writer->SetFileName(fileName.c_str());
        writer->SetInputData(ugrid);
        writer->Write();
        */
        // Add the polygon to a list of polygons
        vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();
        polygons->InsertNextCell(polygon);

        // Create a PolyData
        vtkSmartPointer<vtkPolyData> polygonPolyData = vtkSmartPointer<vtkPolyData>::New();
        polygonPolyData->SetPoints(points);
        polygonPolyData->SetPolys(polygons);
        polygonPolyData->BuildLinks();
       // polygon->Delete();

        // Create a mapper and actor
        //vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
        //mapper->SetInputData(polygonPolyData)
        /*
        vtkNew<vtkExodusIIWriter> exodusWriter;
        exodusWriter->SetFileName("output_ExodusII.e");
       //exodusWriter->SetInputConnection(ugrid->GetOutputPort());
        exodusWriter->SetInputDataObject(ugrid);
        exodusWriter->WriteAllTimeStepsOn();
        exodusWriter->Write();
        */
        vtkNew<vtkDataSetMapper> ugridMapper;
        ugridMapper->SetInputData(ugrid);

        string exodusii_name = "fine_grid_data.e";
        //char p3[20] = "fault.e";
        //vtkSmartPointer<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkNew<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkExodusIIWriter->SetFileName(exodusii_name.c_str());
        vtkExodusIIWriter->WriteOutBlockIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalNodeIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalElementIdArrayOn();
        vtkExodusIIWriter->WriteAllTimeStepsOn();
        //print "#time", vtkTemporalDataSet.GetNumberOfTimeSteps();
        //vtkExodusIIWriter->SetInputData(polygonPolyData);
        vtkExodusIIWriter->SetInputData(ugrid);
        //vtkExodusIIWriter->SetInputConnection(ugridMapper->GetOutputPort());
        vtkExodusIIWriter->Update();
        vtkExodusIIWriter->Write();
        //vtkExodusIIWriter->Delete();

        vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
        actor->SetMapper(ugridMapper);
        actor->GetProperty()->SetColor(colors->GetColor3d("Silver").GetData());

        // Visualize
        vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
        vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
        renderWindow->SetWindowName("Polygon");
        renderWindow->AddRenderer(renderer);
        vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
        renderWindowInteractor->SetRenderWindow(renderWindow);

        renderer->AddActor(actor);
        renderer->SetBackground(colors->GetColor3d("DarkGray").GetData());
        //renderer->SetBackground(0,0,0);
        renderWindow->Render();
        renderWindowInteractor->Start();
    }

    void buildFaultUnstructuredGridB(const string& path, const string& file_name, const string& output_path)
    {

        string file_path = path + file_name;
        fstream fin;
        fin.open(file_path, ios::in);

        if (!fin.is_open())
        {
            cout << "File not found !" << file_path << endl;
        }
        string line;
        double x, y, z;
        int group_id;
        vector<GroupPoint3D> points_v;

        getline(fin, line);
        //while (getline(fin, line))
        while (fin >> x >> y >> group_id >> z) {
            GroupPoint3D gp3d = GroupPoint3D(x, y, z, group_id);
            points_v.push_back(gp3d);
        }
        fin.close();

        ///////////////////////////////////////////////////////////////////////////////
        vtkSmartPointer<vtkNamedColors> colors = vtkSmartPointer<vtkNamedColors>::New();
        vtkNew<vtkUnstructuredGrid> ugrid = vtkNew<vtkUnstructuredGrid>::vtkNew();

        // Setup four points
        vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
        int i = 0;
        for (auto a : points_v) {
            //points->InsertNextPoint(a.xyz);
            points->InsertPoint(i, a.xyz);
            i++;
        }
        // Create the polygon
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        polygon->GetPointIds()->SetNumberOfIds(9); //make a quad

        int index = 0;
        for (int i = 1; i < 10; i++) {
            polygon->GetPointIds()->SetId(index, i);
            polygon->GetPointIds()->SetId(index + 1, i);
            polygon->GetPointIds()->SetId(index + 2, i);

            index = index + 3;
            cout << "index :" << index << " i :" << i << endl;
        }
        for (int i = 0; i < points_v.size(); i = i + 3)
        {
            vtkSmartPointer <vtkIdList > idList_v = vtkSmartPointer<vtkIdList>::New();
            idList_v.Get()->InsertNextId(i);
            idList_v.Get()->InsertNextId(i + 1);
            idList_v.Get()->InsertNextId(i + 2);

            ugrid->InsertNextCell(VTK_LINE, idList_v);
        }
        ugrid->SetPoints(points);
        // Add the polygon to a list of polygons
        vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();
        polygons->InsertNextCell(polygon);

        // Create a PolyData
        vtkSmartPointer<vtkPolyData> polygonPolyData = vtkSmartPointer<vtkPolyData>::New();
        polygonPolyData->SetPoints(points);
        polygonPolyData->SetPolys(polygons);
        polygonPolyData->BuildLinks();

        vtkNew<vtkDataSetMapper> ugridMapper;

        ugridMapper->SetInputData(ugrid);

        
        string exodusii_name = output_path + file_name.substr(0, file_name.size() - 4) + ".e";
        
        //vtkSmartPointer<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkNew<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkExodusIIWriter->SetFileName(exodusii_name.c_str());
        vtkExodusIIWriter->WriteOutBlockIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalNodeIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalElementIdArrayOn();
        vtkExodusIIWriter->WriteAllTimeStepsOn();
       
        //vtkExodusIIWriter->SetInputData(polygonPolyData);
        vtkExodusIIWriter->SetInputData(ugrid);
        //vtkExodusIIWriter->SetInputConnection(ugridMapper->GetOutputPort());
        vtkExodusIIWriter->Update();
        vtkExodusIIWriter->Write();
        vtkExodusIIWriter->Delete();

        vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
        actor->SetMapper(ugridMapper);
        actor->GetProperty()->SetColor(colors->GetColor3d("Silver").GetData());

        // Visualize
        vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
        vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
        renderWindow->SetWindowName("Polygon");
        renderWindow->AddRenderer(renderer);
        vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
        renderWindowInteractor->SetRenderWindow(renderWindow);

        renderer->AddActor(actor);
        renderer->SetBackground(colors->GetColor3d("Salmon").GetData());
        renderWindow->Render();
        renderWindowInteractor->Start();

    }

    void buildQuadWithUnstructuredFineDataGridB(const string& path, const string& file_name, string output_path)
    {
        string file_path = path + file_name;
        fstream fin;
        fin.open(file_path, ios::in);

        if (!fin.is_open())
        {
            cout << "File not found !" << file_path << endl;
        }
        string line;
        double x, y, z;
        int group_id;
        vector<GroupPoint3D> points_v;

        while (fin >> x >> y >> z) {

            GroupPoint3D gp3d = GroupPoint3D(x, y, z);
            points_v.push_back(gp3d);
        }
        fin.close();

        ///////////////////////////////////////////////////////////////////////////////
        vtkSmartPointer<vtkNamedColors> colors = vtkSmartPointer<vtkNamedColors>::New();

        vtkNew<vtkUnstructuredGrid> ugrid = vtkNew<vtkUnstructuredGrid>::vtkNew();
        //vtkNew<vtkUnstructuredGrid> ugrid;
        // Setup four points
        vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
        int i = 0;
        for (auto a : points_v) {
            //cout << a.group_id << " " << x << " " << y << " " << z << endl;
            //points->InsertNextPoint(a.xyz);
            points->InsertPoint(i, a.xyz);
            i++;
        }

        // Create the polygon
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        //vtkNew< vtkPolygon> polygon;
        polygon->GetPointIds()->SetNumberOfIds(46816); //make a quad

        double pts[11704][4];

        int in = 0;
        for (int i = 0; i < 11704; i++) {
            for (int s = 0; s < 4; s++) {
                pts[i][s] = in + s;
            }
            in = in + 4;
        }
        for (int i = 0; i < 11704; i++) {
            for (int j = 0; j < 4; j++) {
                polygon->GetPointIds()->SetId(pts[i][j], i);
            }
        }
        for (int i = 0; i < 11704; i++) {
            vtkSmartPointer <vtkIdList > idList_v = vtkSmartPointer<vtkIdList>::New();
            for (int j = 0; j < 4; j++) {
                idList_v.Get()->InsertNextId(pts[i][j]);
            }
            ugrid->InsertNextCell(VTK_QUAD, idList_v);
        }
        ugrid->SetPoints(points);

       
        // Add the polygon to a list of polygons
        vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();
        polygons->InsertNextCell(polygon);

        // Create a PolyData
        vtkSmartPointer<vtkPolyData> polygonPolyData = vtkSmartPointer<vtkPolyData>::New();
        polygonPolyData->SetPoints(points);
        polygonPolyData->SetPolys(polygons);
        polygonPolyData->BuildLinks();
        //polygon->Delete();

        
        vtkNew<vtkDataSetMapper> ugridMapper;
        ugridMapper->SetInputData(ugrid);

       
        string exodusii_name = output_path + file_name.substr(0, file_name.size() - 4) + ".e";
      
        //vtkSmartPointer<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkNew<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkExodusIIWriter->SetFileName(exodusii_name.c_str());
        vtkExodusIIWriter->WriteOutBlockIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalNodeIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalElementIdArrayOn();
        vtkExodusIIWriter->WriteAllTimeStepsOn();
        
        //vtkExodusIIWriter->SetInputData(polygonPolyData);
        vtkExodusIIWriter->SetInputData(ugrid);
        //vtkExodusIIWriter->SetInputConnection(ugridMapper->GetOutputPort());
        vtkExodusIIWriter->Update();
        vtkExodusIIWriter->Write();
        //vtkExodusIIWriter->Delete();

        vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
        actor->SetMapper(ugridMapper);
        actor->GetProperty()->SetColor(colors->GetColor3d("Silver").GetData());

        // Visualize
        vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
        vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
        renderWindow->SetWindowName("Polygon");
        renderWindow->AddRenderer(renderer);
        vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
        renderWindowInteractor->SetRenderWindow(renderWindow);

        renderer->AddActor(actor);
        renderer->SetBackground(colors->GetColor3d("Salmon").GetData());
        //renderer->SetBackground(0,0,0);
        renderWindow->Render();
        renderWindowInteractor->Start();
    }


    void buildUnstructuredGrid(const string& file_name)
    {

        fstream fin;
        fin.open(file_name, ios::in);

        if (!fin.is_open())
        {
            cout << "File not found !" << file_name << endl;
        }
        string line;
        double x, y, z;
        int group_id;
        vector<GroupPoint3D> points_v;

        getline(fin, line);
        //while (getline(fin, line))
        while (fin >> x >> y >> group_id >> z) {
            GroupPoint3D gp3d = GroupPoint3D(x, y, z, group_id);
            points_v.push_back(gp3d);
        }
        fin.close();

        ///////////////////////////////////////////////////////////////////////////////
        vtkSmartPointer<vtkNamedColors> colors = vtkSmartPointer<vtkNamedColors>::New();
        vtkNew<vtkUnstructuredGrid> ugrid = vtkNew<vtkUnstructuredGrid>::vtkNew();

        // Setup four points
        vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
        int i = 0;
        for (auto a : points_v) {
            //points->InsertNextPoint(a.xyz);
            points->InsertPoint(i, a.xyz);
            i++;
        }

        // Create the polygon
        vtkSmartPointer<vtkPolygon> polygon = vtkSmartPointer<vtkPolygon>::New();
        polygon->GetPointIds()->SetNumberOfIds(22); //make a quad

        int index = 0;
        for (int i = 1; i < 23; i++) {
            polygon->GetPointIds()->SetId(index, i);
            polygon->GetPointIds()->SetId(index + 1, i);
            index = index + 2;
            cout << "index :" << index << " i :" << i << endl;
        }

        for (int i = 0; i < points_v.size(); i = i + 2)
        {
            vtkSmartPointer <vtkIdList > idList_v = vtkSmartPointer<vtkIdList>::New();
            idList_v.Get()->InsertNextId(i);
            idList_v.Get()->InsertNextId(i + 1);
            ugrid->InsertNextCell(VTK_LINE, idList_v);
        }
        ugrid->SetPoints(points);
      
        /*
        string fileName = "result/xmlunstructuredgrid.vtu";
        vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
        //writer.SetFileName('cosipy.vtu')
        writer->SetFileName(fileName.c_str());
        writer->SetInputData(ugrid);
        writer->Write();
        */

        // Add the polygon to a list of polygons
        vtkSmartPointer<vtkCellArray> polygons = vtkSmartPointer<vtkCellArray>::New();
        polygons->InsertNextCell(polygon);

        // Create a PolyData
        vtkSmartPointer<vtkPolyData> polygonPolyData = vtkSmartPointer<vtkPolyData>::New();
        polygonPolyData->SetPoints(points);
        polygonPolyData->SetPolys(polygons);
        polygonPolyData->BuildLinks();

        // Create a mapper and actor
        //vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
        //mapper->SetInputData(polygonPolyData)
        /*
        vtkNew<vtkExodusIIWriter> exodusWriter;
        exodusWriter->SetFileName("output_ExodusII.e");
       // exodusWriter->SetInputConnection(ugrid->GetOutputPort());
        exodusWriter->SetInputDataObject(ugrid);
        exodusWriter->WriteAllTimeStepsOn();
        exodusWriter->Write();
        */
        vtkNew<vtkDataSetMapper> ugridMapper;

        ugridMapper->SetInputData(ugrid);
        
        string exodusii_name = "fault.e";
        //char p3[20] = "fault.e";
        //vtkSmartPointer<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkNew<vtkExodusIIWriter> vtkExodusIIWriter;
        vtkExodusIIWriter->SetFileName(exodusii_name.c_str());
        vtkExodusIIWriter->WriteOutBlockIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalNodeIdArrayOn();
        vtkExodusIIWriter->WriteOutGlobalElementIdArrayOn();
        vtkExodusIIWriter->WriteAllTimeStepsOn();
        //print "#time", vtkTemporalDataSet.GetNumberOfTimeSteps();
        //vtkExodusIIWriter->SetInputData(polygonPolyData);
        vtkExodusIIWriter->SetInputData(ugrid);
        //vtkExodusIIWriter->SetInputConnection(ugridMapper->GetOutputPort());
        vtkExodusIIWriter->Update();
        vtkExodusIIWriter->Write();
        vtkExodusIIWriter->Delete();


        vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
        actor->SetMapper(ugridMapper);
        actor->GetProperty()->SetColor(colors->GetColor3d("Silver").GetData());

        // Visualize
        vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
        vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
        renderWindow->SetWindowName("Polygon");
        renderWindow->AddRenderer(renderer);
        vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
        renderWindowInteractor->SetRenderWindow(renderWindow);

        renderer->AddActor(actor);
        renderer->SetBackground(colors->GetColor3d("Salmon").GetData());
        renderWindow->Render();
        renderWindowInteractor->Start();
    }

    void readExodusFile(string& filename)
    {
        vtkNew<vtkNamedColors> colors;
        // Read Exodus Data
        vtkNew<vtkExodusIIReader> reader;
        reader->SetFileName(filename.c_str());
        reader->UpdateInformation();
        reader->SetTimeStep(10);
        reader->SetAllArrayStatus(vtkExodusIIReader::NODAL,1); // enables all NODAL variables

        // Create Geometry
        vtkNew<vtkCompositeDataGeometryFilter> geometry;
        geometry->SetInputConnection(0, reader->GetOutputPort(0));

        // Mapper
        vtkNew<vtkPolyDataMapper> mapper;
        mapper->SetInputConnection(geometry->GetOutputPort());
        //mapper->SelectColorArray(argv[2]);
        mapper->SetScalarModeToUsePointFieldData();
        mapper->InterpolateScalarsBeforeMappingOn();

        // Actor
        vtkNew<vtkActor> actor;
        actor->SetMapper(mapper);

        // Renderer
        vtkNew<vtkRenderer> renderer;
        renderer->AddViewProp(actor);
        renderer->SetBackground(colors->GetColor3d("DimGray").GetData());

        renderer->GetActiveCamera()->SetPosition(9.0, 9.0, 7.0);
        renderer->GetActiveCamera()->SetFocalPoint(0, 0, 0);
        renderer->GetActiveCamera()->SetViewUp(0.2, -0.7, 0.7);
        renderer->GetActiveCamera()->SetDistance(14.5);

        // Window and Interactor
        vtkNew<vtkRenderWindow> window;
        window->AddRenderer(renderer);
        window->SetSize(600, 600);
        window->SetWindowName("ReadExodusData");

        vtkNew<vtkRenderWindowInteractor> interactor;
        interactor->SetRenderWindow(window);
        interactor->Initialize();

        // Show the result
        window->Render();
        interactor->Start();
    }

    void writeExodusFile()
    {
        vtkNew<vtkTimeSourceExample> timeSource;

        vtkNew<vtkExodusIIWriter> exodusWriter;
        exodusWriter->SetFileName("output.exii");
        exodusWriter->SetInputConnection(timeSource->GetOutputPort());
       // exodusWriter->SetInputDataObject()
        exodusWriter->WriteAllTimeStepsOn();     
       
        exodusWriter->Write();

    }

    vtkPolyData* custom_reader(std::istream& infile) {
        vtkIdType number_of_points, number_of_triangles;
        infile >> number_of_points >> number_of_triangles;
        vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
        points->SetNumberOfPoints(number_of_points);


        for (vtkIdType i = 0; i < number_of_points; i++) {
            double x, y, z;
            infile >> x >> y >> z;
            points->SetPoint(i, x, y, z);

        }
        vtkSmartPointer<vtkCellArray> polys  = vtkSmartPointer<vtkCellArray>::New();
        for (vtkIdType i = 0; i < number_of_triangles; i++)
        {
            vtkIdType a, b, c;
            infile >> a >> b >> c;
            polys->InsertNextCell(3);
            polys->InsertCellPoint(a);
            polys->InsertCellPoint(b);
            polys->InsertCellPoint(c);
        }
        vtkPolyData* polyData = vtkPolyData::New();
        polyData->SetPoints(points);
        polyData->SetPolys(polys);
        return polyData;


    }

    void interploate(const string& file0, const string& file1, const string& prop_file, const string& yougmod) {
        

        std::string variable = "u";
        double range[2] = { 0, 10 };
        ///////////////////////////////////////////////////////////////////////////
        // FILE 0: COARSE MESH WITH SOLUTION 0
        vtkNew<vtkExodusIIReader> coarseReader;
        coarseReader->SetFileName(file0.c_str());
        coarseReader->UpdateInformation();
        coarseReader->SetTimeStep(0);
        coarseReader->SetAllArrayStatus(vtkExodusIIReader::NODAL, 1);
        coarseReader->Update();

        vtkNew<vtkCompositeDataGeometryFilter> coarseGeometry;
        coarseGeometry->SetInputConnection(0, coarseReader->GetOutputPort(0));
        coarseGeometry->Update();

        vtkNew<vtkPolyDataMapper> coarseMapper;
        coarseMapper->SetInputConnection(coarseGeometry->GetOutputPort());
        coarseMapper->SelectColorArray(variable.c_str());
        coarseMapper->SetScalarModeToUsePointFieldData();
        coarseMapper->InterpolateScalarsBeforeMappingOn();
        coarseMapper->SetScalarRange(range);


        //create scalarbar
       /** vtkSmartPointer<vtkLookupTable> lut = vtkSmartPointer<vtkLookupTable>::New();
        lut->SetNumberOfColors(256);
        lut->SetHueRange(0.667, 0);
        //lut->SetRange(range[0], range[1]);  //value below 0.01
        lut->SetRange(0, 1);  //value below 0.01
        lut->Build();
        coarseMapper->SetLookupTable(lut);*/

        vtkNew<vtkActor> coarseActor;
        coarseActor->SetMapper(coarseMapper);
        coarseActor->GetProperty()->SetEdgeVisibility(true);

        vtkNew<vtkRenderer> coarseRenderer;
        coarseRenderer->AddViewProp(coarseActor);

        ///////////////////////////////////////////////////////////////////////////
        // YOUNGMOD 0: YOUNGMOD WITH SOLUTION 0

        vtkNew<vtkExodusIIReader> youngPropReader;

        youngPropReader->SetFileName(yougmod.c_str());
        youngPropReader->UpdateInformation();
        youngPropReader->SetTimeStep(0);
        youngPropReader->SetAllArrayStatus(vtkExodusIIReader::NODAL, 1);
        youngPropReader->Update();

        vtkNew<vtkCompositeDataGeometryFilter> youngGeometry;
        youngGeometry->SetInputConnection(0, youngPropReader->GetOutputPort(0));
        youngGeometry->Update();

        ///////////////////////////////////////////////////////////////////////////
       // PISSON 0: PISSON WITH SOLUTION 0
        vtkNew<vtkExodusIIReader> pissonPropReader;

        pissonPropReader->SetFileName(prop_file.c_str());
        pissonPropReader->UpdateInformation();
        pissonPropReader->SetTimeStep(0);
        pissonPropReader->SetAllArrayStatus(vtkExodusIIReader::NODAL, 1);
        pissonPropReader->Update();

        vtkNew<vtkCompositeDataGeometryFilter> pissonGeometry;
        pissonGeometry->SetInputConnection(0, pissonPropReader->GetOutputPort(0));
        pissonGeometry->Update();

        vtkNew<vtkPolyDataMapper> pissonMapper;
        pissonMapper->SetInputConnection(pissonGeometry->GetOutputPort());
        pissonMapper->SelectColorArray(variable.c_str());
        pissonMapper->SetScalarModeToUsePointFieldData();
        pissonMapper->InterpolateScalarsBeforeMappingOn();
        pissonMapper->SetScalarRange(range);

        vtkNew<vtkActor> pissonActor;
        pissonActor->SetMapper(pissonMapper);
        pissonActor->GetProperty()->SetEdgeVisibility(true);

        vtkNew<vtkRenderer> pissonRenderer;
        pissonRenderer->AddViewProp(pissonActor);

        /////////////////////////////////////////////////////////////////////
        // BASE GRID : 
        vtkNew<vtkExodusIIReader> baseGridReader;

        baseGridReader->SetFileName(file0.c_str());
        baseGridReader->UpdateInformation();
        baseGridReader->SetTimeStep(0);
        baseGridReader->SetAllArrayStatus(vtkExodusIIReader::NODAL, 1);
        baseGridReader->Update();

        vtkNew<vtkCompositeDataGeometryFilter> baseGridGeometry;
        baseGridGeometry->SetInputConnection(0, baseGridReader->GetOutputPort(0));
        baseGridGeometry->Update();

        vtkNew<vtkPolyDataMapper> baseGridMapper;
        baseGridMapper->SetInputConnection(baseGridGeometry->GetOutputPort());
        baseGridMapper->SelectColorArray(variable.c_str());
        baseGridMapper->SetScalarModeToUsePointFieldData();
        baseGridMapper->InterpolateScalarsBeforeMappingOn();
        baseGridMapper->SetScalarRange(range);

        vtkNew<vtkActor> baseGridActor;
        baseGridActor->SetMapper(baseGridMapper);
        baseGridActor->GetProperty()->SetEdgeVisibility(true);

        vtkNew<vtkRenderer> baseGridRenderer;
        baseGridRenderer->AddViewProp(baseGridActor);


        /////////////////////////////////////////////////////////////////////
        // FILE 1: FINE MESH WITH SOLUTION 1

        vtkNew<vtkExodusIIReader> fineReader;
        fineReader->SetFileName(file1.c_str());
        fineReader->UpdateInformation();
        fineReader->SetTimeStep(0);
        fineReader->SetAllArrayStatus(vtkExodusIIReader::NODAL, 1);
        fineReader->Update();

        vtkNew<vtkCompositeDataGeometryFilter> fineGeometry;
        fineGeometry->SetInputConnection(0, fineReader->GetOutputPort(0));
        fineGeometry->Update();
        fineGeometry->GetOutput()->GetPointData()->SetActiveScalars(variable.c_str());

        vtkNew<vtkPolyDataMapper> fineGeometryMapper;
        fineGeometryMapper->SetInputConnection(fineGeometry->GetOutputPort());
        fineGeometryMapper->SelectColorArray(variable.c_str());
        fineGeometryMapper->SetScalarModeToUsePointFieldData();
        fineGeometryMapper->InterpolateScalarsBeforeMappingOn();
        fineGeometryMapper->SetScalarRange(range);

        vtkNew<vtkActor> fineGeometryActor;
        fineGeometryActor->SetMapper(fineGeometryMapper);
        fineGeometryActor->GetProperty()->SetEdgeVisibility(true);

        vtkNew<vtkRenderer> fineGeometryRenderer;
        fineGeometryRenderer->AddViewProp(fineGeometryActor);

        //////////////////////////////////////////////////////////////////////
        // PROJECT SOLUTION FROM FILE 0 to GRID FROM FILE 1

        // Build the structure to interpolate onto
        vtkSmartPointer<vtkUnstructuredGrid>
            coarseInterpolatedGrid = // output to be interpolated on to
            vtkSmartPointer<vtkUnstructuredGrid>::New();

        vtkMultiBlockDataSet* coarseMultiBlock = dynamic_cast<vtkMultiBlockDataSet*>(
            coarseReader->GetOutput()->GetBlock(0));
        coarseInterpolatedGrid->DeepCopy(
            dynamic_cast<vtkUnstructuredGrid*>(coarseMultiBlock->GetBlock(0)));

        vtkNew<vtkStaticPointLocator> locator;
        locator->SetDataSet(fineGeometry->GetOutput());        
        locator->BuildLocator();

        vtkNew<vtkGaussianKernel> kernel;
        kernel->SetSharpness(4);
        kernel->SetKernelFootprintToNClosest();
        kernel->SetNumberOfPoints(10);
        kernel->SetSharpness(4.0);

        // Probe the fine geometry with the course geometry.
        // NOTE: The point arrays will not be interpolated unless PassPointArrays is
        // off.
        vtkNew<vtkPointInterpolator> coarseInterpolator;
        //coarseInterpolator->SetSourceData(fineGeometry->GetOutput()); // Pc data set to be probed by input points P
        coarseInterpolator->SetSourceData(baseGridGeometry->GetOutput());
        coarseInterpolator->SetInputData(coarseGeometry->GetOutput());

        coarseInterpolator->SetKernel(kernel);
        coarseInterpolator->SetLocator(locator);
        coarseInterpolator->SetNullPointsStrategyToClosestPoint();
        coarseInterpolator->PassPointArraysOff();
        coarseInterpolator->Update();

        string filename = "inteplator_grid_01.vtk";
        //writer = vtk.vtkXMLImageDataWriter()
        //writer.SetFileName(filename)
        //writer.SetInputData(interpolator.GetOutput())
        //writer.Write()
        //pcl::PolygonMesh mesh;
        vtkSmartPointer <vtkPolyDataWriter> polyWriter = vtkSmartPointer<vtkPolyDataWriter>::New();
        //vtkNew<vtkXMLPolyDataWriter> polyWriter;
        polyWriter->SetFileName(filename.c_str());
        polyWriter->SetInputData(coarseInterpolator->GetOutput());
        polyWriter->Write();
       
        vtkNew<vtkDataSetMapper> coarseInterpolatorMapper;
        coarseInterpolatorMapper->SelectColorArray(variable.c_str());
        coarseInterpolatorMapper->SetScalarModeToUsePointFieldData();
        coarseInterpolatorMapper->InterpolateScalarsBeforeMappingOn();
        coarseInterpolatorMapper->SetScalarRange(range);

        vtkNew<vtkActor> coarseInterpolatorActor;
        coarseInterpolatorActor->SetMapper(coarseInterpolatorMapper);
        coarseInterpolatorActor->GetProperty()->SetEdgeVisibility(true);

        vtkNew<vtkRenderer> coarseInterpolatorRenderer;
        coarseInterpolatorRenderer->AddActor(coarseInterpolatorActor);
     
        // Set the active scalar for the two inputs.
        // NOTE: InterpolateDataSetAttibutes does not interpolate field data.
        // To interpolate field data it must be added as the active scalar
        coarseInterpolatedGrid->GetPointData()->SetActiveScalars(variable.c_str());
        coarseInterpolator->GetOutput()->GetPointData()->SetActiveScalars( variable.c_str());
        vtkNew<vtkInterpolateDataSetAttributes> coarseInterpolateAttributes;
        coarseInterpolateAttributes->AddInputData(0, coarseInterpolatedGrid);
        coarseInterpolateAttributes->AddInputData(0, coarseInterpolator->GetOutput());
        coarseInterpolateAttributes->AddInputConnection(baseGridGeometry->GetOutputPort());
        coarseInterpolateAttributes->AddInputData(0,pissonGeometry->GetOutputDataObject(0));
        coarseInterpolateAttributes->AddInputData(0, youngGeometry->GetOutputDataObject(0));
       
        coarseInterpolateAttributes->SetT(0.5);
        coarseInterpolateAttributes->Update();

        vtkNew<vtkDataSetMapper> coarseInterpolateAttibutesMapper;
        coarseInterpolateAttibutesMapper->SetInputConnection( coarseInterpolateAttributes->GetOutputPort());
        coarseInterpolateAttibutesMapper->SelectColorArray(variable.c_str());
        coarseInterpolateAttibutesMapper->SetScalarModeToUsePointFieldData();
        coarseInterpolateAttibutesMapper->InterpolateScalarsBeforeMappingOn();
        coarseInterpolateAttibutesMapper->SetScalarRange(range);

        vtkNew<vtkActor> coarseInterpolateAttributesActor;
        coarseInterpolateAttributesActor->SetMapper(coarseInterpolateAttibutesMapper);
        coarseInterpolateAttributesActor->GetProperty()->SetEdgeVisibility(true);


        vtkNew<vtkRenderer> coarseInterpolateAttributesRenderer;
        coarseInterpolateAttributesRenderer->AddActor(coarseInterpolateAttributesActor);

        /////////////////////////
        //  Window and Interactor

        vtkNew<vtkRenderWindow> window;

        // Use the same camera for each renderer
        fineGeometryRenderer->SetActiveCamera(coarseRenderer->GetActiveCamera());
        coarseInterpolatorRenderer->SetActiveCamera(
            coarseRenderer->GetActiveCamera());
        coarseInterpolateAttributesRenderer->SetActiveCamera(
            coarseRenderer->GetActiveCamera());

        // Set each renderer color
        vtkNew<vtkNamedColors> colors;
       // coarseRenderer->SetBackground(colors->GetColor3d("Gainsboro").GetData());
        //coarseInterpolatorRenderer->SetBackground(
        //    colors->GetColor3d("LightGrey").GetData());
        coarseInterpolateAttributesRenderer->SetBackground(
            colors->GetColor3d("Silver").GetData());
        fineGeometryRenderer->SetBackground(colors->GetColor3d("DarkGray").GetData());

        // Set the viewport for each renderer
        coarseInterpolateAttributesRenderer->SetViewport(0, 0, 0.5, 1);
        fineGeometryRenderer->SetViewport(0.5, 0, 1, 1);
        //coarseInterpolateAttributesRenderer->SetViewport(0.5, 0, 0.75, 1);
        //fineGeometryRenderer->SetViewport(0.75, 0, 1, 1);

        //Add the renderers
        //window->AddRenderer(coarseRenderer);             // low res
        //window->AddRenderer(coarseInterpolatorRenderer); // interpolated geometry
        window->AddRenderer(
            coarseInterpolateAttributesRenderer);  // interpolated attributes
        window->AddRenderer(fineGeometryRenderer); // hi res with point data

        // The size the window to hold 4 viewports of width 320
        window->SetSize(1280, 320);
        window->SetWindowName("InterpolateFieldDataDemo");

        // Create the interactor
        vtkNew<vtkRenderWindowInteractor> interactor;
        interactor->SetRenderWindow(window);
        window->Render();
        interactor->Initialize();

        // Reset the camera to a default view, shared by all
        coarseRenderer->ResetCamera();
        interactor->Start();
        ///////////////////////////////////////////////////////

       //pcl::PolygonMesh mesh;
       //vtkSmartPointer <vtkPolyDataWriter> polyWriter = vtkSmartPointer<vtkPolyDataWriter>::New();
        /*vtkNew<vtkXMLPolyDataWriter> polyWriter;
        polyWriter->SetFileName("inteplator_grid.vtk");
        vtkSmartPointer<vtkPolyData> polyInputData;
        polyInputData = coarseInterpolateAttributes->GetPolyDataOutput();
        polyWriter->SetInputData(polyInputData);
        
        //polyWriter->SetInputData(coarseInterpolateAttributes->GetOutput()->GetPointData());

        polyWriter->Write();*/

        ///////////////////////////////////////////////////////
    }

    void interploateLetter() {
        double t = 2.0;
       
        // make the letter v
        vtkNew<vtkVectorText> letterV;
        letterV->SetText("v");

        // read the geometry file containing the letter t
        vtkNew<vtkVectorText> letterT;
        letterT->SetText("t");

        // read the geometry file containing the letter k
        vtkNew<vtkVectorText> letterK;
        letterK->SetText("k");

        // create implicit models of each letter
        vtkNew<vtkImplicitModeller> blobbyV;
        //letterV->GetOutputPort
        blobbyV->SetInputConnection(letterV->GetOutputPort());
        blobbyV->SetMaximumDistance(.2);
        blobbyV->SetSampleDimensions(50, 50, 12);
        blobbyV->SetModelBounds(-0.5, 1.5, -0.5, 1.5, -0.5, 0.5);

        vtkNew<vtkImplicitModeller> blobbyT;
        blobbyT->SetInputConnection(letterT->GetOutputPort());
        blobbyT->SetMaximumDistance(.2);
        blobbyT->SetSampleDimensions(50, 50, 12);
        blobbyT->SetModelBounds(-0.5, 1.5, -0.5, 1.5, -0.5, 0.5);

        vtkNew<vtkImplicitModeller> blobbyK;
        blobbyK->SetInputConnection(letterK->GetOutputPort());
        blobbyK->SetMaximumDistance(.2);
        blobbyK->SetSampleDimensions(50, 50, 12);
        blobbyK->SetModelBounds(-0.5, 1.5, -0.5, 1.5, -0.5, 0.5);

        // Interpolate the data
        vtkNew<vtkInterpolateDataSetAttributes> interpolate;
        interpolate->AddInputConnection(blobbyV->GetOutputPort());
        interpolate->AddInputConnection(blobbyT->GetOutputPort());
        interpolate->AddInputConnection(blobbyK->GetOutputPort());
        interpolate->SetT(0.0);

        // extract an iso surface
        vtkNew<vtkContourFilter> blobbyIso;
        blobbyIso->SetInputConnection(interpolate->GetOutputPort());
        blobbyIso->SetValue(0, 0.1);

        // map to rendering primitives
        vtkNew<vtkPolyDataMapper> blobbyMapper;
        blobbyMapper->SetInputConnection(blobbyIso->GetOutputPort());
        blobbyMapper->ScalarVisibilityOff();

        vtkNew<vtkNamedColors> colors;

        // now an actor
        vtkNew<vtkActor> blobby;
        blobby->SetMapper(blobbyMapper);
        blobby->GetProperty()->SetDiffuseColor(
            colors->GetColor3d("Banana").GetData());
        blobby->GetProperty()->EdgeVisibilityOn();

        // Create the RenderWindow, Renderer and both Actors
        //
        vtkNew<vtkRenderer> renderer;
        vtkNew<vtkRenderWindow> renderWindow;
        renderWindow->AddRenderer(renderer);
        vtkNew<vtkRenderWindowInteractor> interactor;
        interactor->SetRenderWindow(renderWindow);

        vtkNew<vtkCamera> camera;
        camera->SetClippingRange(0.265, 13.2);
        camera->SetFocalPoint(0.539, 0.47464, 0);
        camera->SetPosition(0.539, 0.474674, 2.644);
        camera->SetViewUp(0, 1, 0);
        camera->Azimuth(30.);
        camera->Elevation(30.);
        renderer->SetActiveCamera(camera);

        //  now  make a renderer and tell it about lights and actors
        renderWindow->SetSize(640, 480);
        renderWindow->SetWindowName("Morph3D");

        renderer->AddActor(blobby);
        renderer->SetBackground(colors->GetColor3d("SlateGray").GetData());
        renderWindow->Render();
        for (unsigned i = 0; i < 100; ++i)
        {
            interpolate->SetT(i / 99.0 * 2.0);
            interpolate->Modified();
            renderWindow->Render();
            vtksys::SystemTools::Delay(100);
        }
        interpolate->SetT(t);
        renderWindow->Render();

        // interact
        interactor->Start();

    }
};

